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Advancements in the synthesis of faceted nanoparticles and colloids have spurred interest in the 
phase behavior of polyhedral shapes. Regular tetrahedra have attracted particular attention because 
they prefer local symmetries that are incompatible with periodicity. Two dense phases of regular 
tetrahedra have been reported recently. The densest known tetrahedron packing is achieved in a 
crystal of triangular bipyramids (dimers) with packing density 4000/4671 ~ 85.63%. In simulation 
a dodecagonal quasicrystal is observed; its approximant, with periodic tiling (3.4.3 2 .4), can be 
compressed to a packing fraction of 85.03%. Here, we show that the quasicrystal approximant 
is more stable than the dimer crystal for packing densities below 84% using Monte Carlo computer 
simulations and free energy calculations. To carry out the free energy calculations, we use a variation 
of the Frenkel-Ladd method for anisotropic shapes and thermodynamic integration. The enhanced 
stability of the approximant can be attributed to a network substructure, which maximizes the free 
volume (and hence the wiggle room) available to the particles and facilitates correlated motion of 
particles, which further contributes to entropy and leads to diffusion for packing densities below 
65%. The existence of a solid-solid transition between structurally distinct phases not related by 
symmetry breaking - the approximant and the dimer crystal- is unusual for hard particle systems. 



I. INTRODUCTION 

The self-assembly of nanoparticles into ordered struc- 
tures is governed by interaction and shape anisotropy [1]. 
Anisotropic particles are capable of stabilizing complex 
phases by entropy alone. Such structures can have poten- 
tially interesting optical and electrical properties yet to 
be fully investigated [2-6]. Among anisotropic particles, 
tetrahedra are promising for assembling unusual struc- 
tures because of their simplicity as well as their lack of 
inversion symmetry. When arranged face-to-face, tetra- 
hedra form configurations with five-fold or icosahedral 
symmetries that are incompatible with periodicity. This 
results in geometric frustration and renders the assem- 
bly of tetrahedra more challenging than assembling other 
shapes. Various types of nano-tetrahedra have recently 
been synthesized from noble metals [7, 8] and crystalline 
silicon [9, 10]. Micron-size colloidal tetrahedra made of 
colloidal spheres have also been reported [11]. In certain 
cases, these tetrahedra may be treated as hard particles. 

Particles whose interactions are dominated by repul- 
sion can be modeled to first approximation as hard parti- 
cles. Since all permissible configurations of such systems 
are of identical energy, entropic effects govern their phase 
behavior. Classic examples of entropy- driven phase tran- 
sitions are the isotropic-to-nematic transition for hard 
thin rods [12] and the crystallization of hard spheres into 
close-packed structures upon compression [13]. Entropy 
drives these particles to order, because doing so will in- 
crease the number of configurations accessible to the sys- 
tem. In other words, the increase in macroscopic (visible) 
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order is accompanied by an increase in microscopic disor- 
der (the number of microstates) [14]. The origin of order- 
ing can also be explained by considering the underlying 
thermodynamics of hard particle systems. In the limit of 
infinite pressure, the Gibbs free energy G = PV — ST is 
dominated by the PV term, which means that the dens- 
est packing will be ultimately stable at sufficiently high 
pressures. To date, all known maximally dense packings 
of hard shapes are ordered [15]. 

Although the phase behavior of hard spheres has been 
investigated extensively [16], many fewer studies have 
been done on other hard shapes [17-25]. A key fea- 
ture of the reported phase diagrams is the occurrence 
of symmetry-breaking phase transitions (first and second 
order) in which the symmetry group of the high-density 
phase is a subgroup of the symmetry group of the low- 
density phase (see for example the phase transitions in 
[25]). This means that the compression of the isotropic 
fluid results in an increase of structural complexity by 
breaking at least one symmetry per transformation. For 
instance, hard cubes form a cubatic liquid crystal before 
crystallizing into a simple cubic lattice. In both the liquid 
crystal and the cubic crystal the rotational symmetry is 
broken while the translational symmetry is only broken 
in the crystal and is present in the cubatic phase [23]. 

The problem of assembling and packing hard tetra- 
hedra has drawn significant attention over the last few 
years [26-36] and two competing phases have been re- 
ported in the high-density regime. The densest known 
packing of regular tetrahedra is a structurally simple 
double-triangular bipyramid crystal with packing density 
<j) = 4000/4671 ~ 85.63% obtained from analytical con- 
struction and supported by numerical simulation [32, 34]. 
It is obtained through optimizing an earlier monoclinic 
crystal discovered by Kallus et al [30, 33]. We refer to 
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this structure as the dimer crystal throughout this work 
since the packing is characterized by pairs of tetrahedra 
incorporated into triangular bipyramids (Fig. la). 

Despite its stability in the limit of infinite pressure, 
simulations show that the dimer crystal does not form 
from the fluid except for systems of 16 or fewer parti- 
cles [34]. Instead, a dodecagonal quasicrystal sponta- 
neously assembles at packing densities close to 50% and 
above [29]. Structurally, the quasicrystal is significantly 
more complicated than the dimer phase; tetrahedra are 
arranged into rings that are further capped with pentag- 
onal dipyramids (PDs). The rings and PDs are stacked 
in logs parallel to the ring axis, which in projection form 
the vertices of a planar tiling of squares and triangles 
(Fig. lb). Additional particles- referred to as intertitials- 
appear in the space between the neighboring logs. It is 
noteworthy that the entire structure is a network of in- 
terpenetrating PDs spanning all particles in the system. 
A periodic approximant of the quasicrystal, i.e. a crys- 
tal approximating the structure of the quasicrystal on a 
local level, with the (3.4.3 2 .4) Archimedean tiling and 82 
tetrahedra per unit cell compresses up to = 85.03%, 
only slightly less dense than the dimer crystal [29]. In 
this paper we demonstrate that the approximant is more 
stable than the dimer crystal up to very high pressures 
and that the system prefers the dimer crystal thermody- 
namically only at packing densities exceeding 84%. 

Here we carry out a detailed investigation of the phase 
behavior of hard tetrahedra from the fluid up to the dens- 
est packing. In contrast to previously studied systems of 
hard particles, the phase diagram of tetrahedra entails 
a non-symmetry-breaking solid-solid transition. We con- 
firm its existence by Monte Carlo simulation and free en- 
ergy calculations and discuss the origin of the transition. 
The present study complements previous works on hard 
tetrahedra which studied some aspects of the equation of 
state [29, 37, 38] as well as dense packings [33, 34], and 
extends those to provide a complete picture of the phase 
diagram. By comparing the results of self-assembly sim- 
ulations to those obtained from free energy calculations, 
we assess the likelihood of various candidate phases to be 
observed both in simulations and in experiments of hard 
tetrahedra. 

The paper is organized as follows. The simulation 
methods as well as technical details of the free energy 
and free volume calculations are presented in Section II. 
In Section III A, the thermodynamics of the dimer phase 
is reported. The thermodynamics of the quasicrystal and 
its approximant follows in Section IIIB. The results of 
free energy calculations are presented in Section IIIC. A 
computer experiment in which the dimer crystal spon- 
taneously transforms into the quasicrystal at <\> = 50% 
is reported in Section HID. The origin of the stability 
of the approximant over the dimer crystal at experimen- 
tally realizable densities is discussed in Section III E and 
discussions and concluding remarks are provided in Sec- 
tion IV. 




FIG. 1: Dense packings of tetrahedra proposed by geomet- 
ric construction and computer simulation, (a) The densest 
known tetrahedron packing [34] is a crystal with four particles 
per unit cell forming two triangular bipyramids (or 'dimers') 
shown in green and blue. The unit cell breaks the three-fold 
symmetry of the dimers. (b) In simulation, the particles form 
rings of twelve tetrahedra (red) capped by pentagonal dipyra- 
mids (green) together with interstitial tetrahedra (blue) . The 
rings stack in logs, which arrange to form the vertices of a pla- 
nar square-triangle tiling. Tilings observed in simulation are 
quasiperiodic, but a denser packing is obtained with the pe- 
riodic (3.4.3 2 .4) Archimedean tiling, which is an approximant 
of the quasicrystal [29] . 



II. METHODS 

Simulations of N hard, regular tetrahedra are carried 
out in the isochoric (NVT) ensemble and the isobaric 
(NPT) ensemble using a Monte Carlo algorithm. For- 
bidden overlaps of tetrahedra are determined using the 
separating axis theorem as explained in detail in Ref. [29] . 
N particle trial moves are executed per Monte Carlo cy- 
cle. Each trial move can be a trial translation or a trial 
rotation chosen with equal probability. In the isobaric 
simulations, an additional box trial move is also per- 
formed where the size and shape of the simulation box 
are changed. The edge length of a tetrahedron, a, is cho- 
sen as the unit length of the system. The effective pres- 
sure P* = Pa 3 /ksT is measured in dimensionless units. 
Maximum steps sizes are adjusted occasionally to allow 
for a target acceptance probability of 30% and periodic 
boundary conditions are applied in all three dimensions. 



A. Equation of state 

Equations of state, 0(P*), are calculated with isobaric 
simulations. Changes in the Gibbs free energy within a 
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single phase are obtained via thermodynamic integration: 



C?2 — G\ 
Nk R T 



dp 



(i) 



where Vp = cr 3 y/2/12 is the volume of a tetrahedron. 

Simulations are carried out in the pressure range 50 < 
P* < 4000 for the dimer crystal (4 x 6 x 6 x 6 = 864 
tetrahedra), quasicrystal (8,000 tetrahedra) assembled 
from the fluid and compressed to a packing density up 
to 83.36%, and the approximant (82 x 2 x 2 x 3 = 984 
tetrahedra). 



B. Pressure estimation 

The acceptance probability of trial volume changes 
is an estimator of the pressure in Monte Carlo simula- 
tions [39]. Consider a trial expansion that increases the 
volume from V to V + AV. To fulfill detailed balance, 
the acceptance probability of the volume change is given 
by the Boltzmann factor, 



P B = exp 



P*AV 



N In 1 



AV 



} 
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On the other hand, a trial compression that decreases 
the volume from V to V — AV is accepted if and only 
if no overlap is generated by the trial volume change. 
Let Pno be the probability to generate an overlap in the 
trial compression. For small AV and in equilibrium the 
probabilities are equal, Pno = Pb, and we can solve for 
the pressure: 
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Here, p^o = P^/o' 1S the probability of a single particle 
not having any overlap with any other particle after the 
trial compression that decreases the volume by AV. 



C. Free energy calculations 

1. Frenkel-Ladd method for anisotropic hard particles 

The free energy of a (quasi-) crystal is calculated us- 
ing the Frenkel-Ladd method [39, 40] by transforming it 
reversibly into an Einstein crystal, which serves as a ref- 
erence structure with known free energy. In the Einstein 
crystal, each particle is tethered to its average lattice po- 
sition via harmonic springs. Although originally devel- 
oped for spherical particles, this method can be extended 
to particles with rotational degrees of freedom, such as 
tetrahedra. Additional springs are needed to tether the 
orientations of the particles to their average orientations 
in the lattice. Alternative extensions of the Frenkel-Ladd 
method to systems of particles with rotational degrees of 
freedom can be found in the literature [41]. 



We describe the configuration of a tetrahedron by 
(r, q), with r being its center of mass position and q the 
unit quaternion describing its orientation. The potential 
energy of of the corresponding Einstein crystal can then 
be expressed as: 
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where r^o and q^ are the reference position and the ref- 
erence orientation of the i-th particle in the crystal. The 
constant c allows us to adjust the relative strength of 
the rotational springs and does not affect the computed 
free energy differences. All the results in this study are 
obtained using a value of c = 1/2; we tested that us- 
ing other values of c does not affect the outcome of the 
calculations. 

Each system is transformed to the Einstein crystal 
along a reversible path parameterized by 7 G [0, 7 max ] 
using the isochoric-isothermal (NVT) ensemble and the 
Hamiltonian 



(5) 



The hard particle system with Hamiltonian Hhard cor- 
responds to 7 = 0, while in the limit 7 —> 00 the Ein- 
stein crystal is obtained. In practice, we can stop at 
a sufficiently large value of 7 max when the springs are 
strong enough to suppress any particle collisions. The 
Helmholtz free energy difference AA = A^i n — Ahard be- 
tween the reference Einstein crystal and the hard particle 
system is given by: 



AA 



(C/) 7 d 7 (6) 



Note that the Frenkel-Ladd method can only be used 
if there is no translational or rotational diffusion in the 
system; otherwise the ensemble average {U) 1 will not be 
well-defined for small values of 7. 

In our simulations, the system is held for 2 x 10 5 Monte 
Carlo cycles at each 7 value during which {U) 1 is eval- 
uated. The integral in equation (6) is then computed 
numerically. This allows us to determine the Gibbs free 
energy G = A + PV of the dimer (D) and the approxi- 
mant (A) in the range 250 < P* < 600 where no config- 
urational rearrangements are observed. The free energy 
difference AG = Gd — Ga is extrapolated to pressures 
outside this range using thermodynamic integration in 
addition to the Frenkel-Ladd method [42]: 
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2. Fluid-solid transition 



We determine the melting pressure P^ by calculating 
the absolute free energies of the solid and fluid. For suf- 
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ficiently large values of 7, the Helmholtz free energy of 
the Einstein crystal is given by [39]: 
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where A = hj (27rmk B T) ' is the de Broglie wavelength. 
iVgym is the number of quaternions corresponding to ori- 
entations that are symmetry-equivalent, which is twice 
the order of the rotation group of the particle. The 
factor 2 arises from the fact that quaternions are inher- 
ently degenerate in describing the orientation i.e. q and 
— q correspond to the same rotation matrix. For a non- 
symmetric particle, the rotation group will have one ele- 
ment (identity) only and N sym = 2. Here, for tetrahedra, 
the rotation group has twelve elements, so N sym — 24. 
The first and the second terms are configurational con- 
tributions resulting from the translational and rotational 
springs. The last term corresponds to momentum contri- 
butions due to translational degrees of freedom. Momen- 
tum contributions due to rotational degrees of freedom 
are identical for the fluid and the solid and are therefore 
not included here. 

The Gibbs free energy of an ideal gas, which approxi- 
mates a real gas in the limit of infinite dilution, is 
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The free energy of the fluid phase is then obtained from 
thermodynamic integration [43]: 
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We calculate Gfl U id(P*) using the equation of state for a 
system of N = 4,096 tetrahedra for 0.01 < P* < 60. 



ously without overlapping with its neighbors while keep- 
ing all the other particles fixed [44]. The definition gen- 
eralizes to anisotropic particles with rotational degrees of 
freedom where free volume Vf is now the volume of the 
largest subset of configurational space connected to the 
origin that can be accessed by a given particle while fixing 
the positions and orientations of all other particles [45]: 
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Here, J(r, q) is the indicator function of motions (r, q) 
consisting of a translation by r and a rotation by q and 
connected to the origin. / is unity if the particle does not 
overlap with any other particle and zero otherwise. Due 
to the inherent periodicity of rotational motion, the free 
volume of an anisotropic particle has generally a more 
complicated topology compared to the free volume of a 
sphere. Here we calculate free volumes at high densities 
where the free volume is simply connected. 



1. Shooting method 

We calculate the free volume of a particle using a 
method we call the shooting method. Let (u, v) corre- 
spond to a unit vector in the six dimensional configura- 
tion space and suppose that particle i is 'shot' in this 
direction until it hits another particle. The 'shooting 
distance' is the smallest value of a for which the particle 
first overlaps with its neighbors if translated by au and 
oriented according to the quaternion (q i +av)/||q i +av||. 

A lower bound for the free volume can be obtained 
by averaging over a sufficiently large number N s of shots 
with shot distances otj along randomly chosen directions: 
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3. Finite size effects 

To ensure the system sizes we use are free of finite 
size effects, we calculate the Gibbs free energy difference 
between the dimer crystal and the approximant AG = 
Gd — Ga, using equation (8): 
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For the particle numbers used in the free energy cal- 
culations, N D = 864, N A = 984, 7 max = 4 x 10 6 and 
a/A = V2tt, the error in AG is on the order of 10 -3 k B T, 
which is negligible for the present purposes. 



D. Free volume calculations 

The free volume of a hard sphere is the volume of the 
region of space in which the sphere can be moved continu- 



Here 7r 3 /6 is the volume of the six-dimensional unit 
sphere. Note that the periodic topology and the cur- 
vature of the six-dimensional configuration space are ig- 
nored, which is acceptable at high packing densities be- 
cause ||Aq|| < 1. 

Eq. (13) is a lower bound for concave free volumes, 
because shooting only allows access to the parts of the 
free volume connected to the origin by a straight line. 
Non-convex free volumes can arise from sliding collisions 
which, however, become increasingly rare at high packing 
densities. In fact, as we will show now for tetrahedra, 
the shooting method is accurate for high enough packing 
densities. 



2. Binning method 

To estimate the amount of error in the shooting 
method introduced by non- convexity, we using the al- 
ternative binning method which corresponds to a Monte 
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Carlo integration of the free volume. The configuration 
space of a given particle is partitioned into A^ins small 
radial bins of volume V bm . We perform N t random ghost 
trial moves per bin to average out the orientational de- 
grees of freedom and determine the number 7V NO of trial 
moves not leading to an overlap. Free volume can then 
be estimated from: 



Dimer I, Di 



Dimer II, Dn 



Dimer III, Dm 



v f 
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Nt 
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j = l 



NO 



(14) 



Binning is much slower than shooting and might overesti- 
mate the free volume, if a trial move discovers an area of 
configuration space without overlap, but not connected 
to the original particle position. We find that the average 
of the logarithms of the free volumes calculated from the 
shooting method and the binning method agree within a 
relative error of 1(T 2 for all densities 6 > 70%. 



3. Mean-field approximation 

The distribution of free volumes is related to the en- 
tropy of a hard particle system in the mean-field approx- 
imation. If we assume that free volumes of neighbor- 
ing particles are uncorrelated, then the partition func- 
tion of the system is expressed as Q m f = Y\4=\ v f,i an d 
the Helmholtz free energy as A^/NksT = —(IriVf). 
The thermodynamically relevant quantity is therefore the 
mean-log average of free volumes: 



V /,ML 



exp(ln Vf) 



(15) 



which will be used in the rest of this study instead of the 
simple average (vf). 



III. RESULTS 

A. Symmetrization of the dimer packing on 
decreasing pressure 

We construct the dimer crystal analytically [34] and 
slowly expand it by reducing the pressure. The crystal re- 
mains stable during the simulation for pressures P* > 60 
while at lower pressures it melts abruptly. No hysteresis 
is observed in the equation of state (Fig. 2a), if the de- 
compression is stopped before melting and the system is 
re-compressed. This suggests that the system remains at 
least in metastable equilibrium over this range of pres- 
sures and densities. 

The compressibility k = (l/<fi)(d<j)/dP*) (Fig. 2b) re- 
veals a complicated phase behavior, with an anomalous 
peak indicative of a second-order phase transition ap- 
pearing at around P* = 90. We verify that this is a 
displacive phase transition; i.e. it only involves a lattice 
distortion and the particles in the lower density phase 
still remain in dimers. Analyzing the lengths of the vec- 
tors spanning the simulation box and the angles between 
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FIG. 2: Symmetrization of the dimer crystal, (a) The equa- 
tion of state shows no hysteresis between compression and 
expansion, (b) A peak is observed in the compressibility near 
pressure P* = 90 at a second order displacive phase tran- 
sition. The (c) three box angles and (d) box lengths, ob- 
tained by sorting the angles and lengths and then averaging 
the sorted values, are plotted as a function of pressure. We 
observe two transitions, from triclinic (Dm) to monoclinic 
(D/i) to rhombohedral (Dj). The phase Dm is thermody- 
namically stable, Du and Dj are metastable. 



them (Figs. 2c, d) indicates that the transition takes place 
in two stages. While in the lower density phase D/ 
(P* < 90) all lengths and angles are equivalent, they 
are completely split only in the phase Dm (P* > 220). 
There is also an intermediate phase Djj (90 < P* < 220) 
in which only two of the lengths and angles are still de- 
generate. The symmetrization of the lattice therefore fol- 
lows the sequence: triclinic (Dm) — >> monoclinic (D/j) 
— >• rhombohedral (Dj). 

It is known that the three-fold symmetry of the dimers 
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must be broken to achieve optimal bulk packing [33, 34], 
and we observe this in the sequence of transitions. We 
note that Djj was initially reported by Kallus et at. as 
a candidate for the densest packing of tetrahedra [33]. 
Its maximum packing density is only 0.2% lower than 
the maximum packing density of D///, the structure pre- 
dicted by Chen et al. [34]. Note also that the integrated 
area under the peak is a measure of the difference in 
packing densities. This explains the missing peak in the 
compressibility for the transition Dm — »• Djj. In con- 
trast, the difference in maximum packing densities for 
the transition D// — >> Dj is much larger, and of the order 
of a few percent. 



B. 



Comparison of the quasicrystal and its (3.4.3 2 .4) 
approximant 



The equations of state of the quasicrystal, the approxi- 
mant, and the dimer crystal are presented in Fig. 3a. We 
observe that the approximant is not only denser than 
the quasicrystal at all pressures above the melting tran- 
sition, it also melts at lower pressure. These observations 
together with Eq. (1) suggest that the quasicrystal is gen- 
erally less stable than the approximant. 

Further evidence for the stability of the approximant 
over the quasicrystal is obtained through constructing 
higher order approximant s, i.e. approximant s that have 
larger unit cells than the (3.4.3 2 .4) approximant, and 
comparing their equations of state with the quasicrys- 
tal and the approximant. For this purpose, we construct 
the second-order approximant with a unit cell contain- 
ing 1,142 tetrahedra using an inflation operation [46] and 
compute its equation of state near the transition region. 

As observed in Fig. 4, the second approximant is denser 
than the densest quasicrystal that formed in our simula- 
tions but less dense than the first approximant. Neither 
structure is expected to have a significant entropic ad- 
vantage over others since tetrahedra experience similar 
local environments in all these structures. It is there- 
fore safe to conclude that the first approximant is more 
stable than the quasicrystal and the second approximant 
because of its higher density. Higher-order approximants 
can be constructed similarly using inflation symmetry; 
however, such approximants will have very large unit 
cells with tens of thousands of particles. Based on the 
observed trend, we expect higher-order approximants to 
become successively less dense but still denser than the 
quasicrystal. 

The question of comparing the relative thermodynamic 
stability of quasicrystals and their approximants plagues 
nearly all reports of new quasicrystals in the literature. 
The difficulty in obtaining perfect quasicrystals in ex- 
periments and simulations, along with the slow kinetics 
that would be involved in the transformation of even an 
imperfect quasicrystal to any of its approximants, con- 
founds attempts to address quasicrystal stability. In this 
spirit, we remark that the quasicrystal configuration used 




500 1000 1500 2000 2500 3000 3500 4000 
Pressure 



FIG. 3: Thermodynamic stability of the dimer crystal, (a) 
The equation of state for the dimer crystal, the approximant, 
and the quasicrystal shows that the dimer crystal is the dens- 
est packing for P* > 700. The approximant is always denser 
than the quasicrystal. Error bars are smaller than the size of 
the symbols. Insets show the equations of state in the melting 
region as well as near P* = 700 where the dimer crystal first 
becomes denser than the approximant. (b) The Gibbs free 
energy difference between the dimer crystal and the approxi- 
mant AG/Nk B T = (G D ~ G A )/Nk B T calculated using ther- 
modynamic integration and the Frenkel-Ladd method. The 
dimer crystal is stable only at very high pressures. 
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FIG. 4: Equations of state of the quasicrystal and the first 
and second approximants computed from NPT simulations. 
The second approximant (with a unit cell of 1,142 particles) 
is less dense than the first approximant (with a 82-particle 
unit cell). 
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FIG. 5: The critical packing fractions for the approximant to 
dimer transition can be calculated via the common tangent 
construction from the Helmholtz free energies of the approx- 
imant (red) and the dimer crystal (blue). 

in this study is obtained in simulation and an ideal, per- 
fect quasicrystal might be slightly denser. The struc- 
ture of such an ideal quasicrystal, however, is unknown. 
A denser quasicrystal would shift the curve in Fig. 3a 
slightly upwards, and hence make the quasicrystal ther- 
modynamically more stable than the approximant in a 
narrow region close to melting. Based on all evidence 
however, we use the (3.4.3 2 .4) approximant as the most 
stable quasicrystal-like structure for free energy and free 
volume calculation purposes. 

C. Relative thermodynamic stability 

The Gibbs free energy difference between the dimer 
and the approximant is calculated using the method de- 
scribed in Sec. II C 1. We find that the dimer crystal is 
stable only for pressures above P* = 3780 ±60 (Fig. 3b), 
while the approximant is favored below P*. At the criti- 
cal pressure, the approximant and the dimer crystal have 
packing densities of (84.0 ±0.1)% and (84.6 ±0.1)% 
respectively. The transition densities can be alterna- 
tively calculated from the Helmholtz energy using the 
common-tangent construction (Fig. 5). P* is signifi- 
cantly higher than the melting pressure for the approx- 
imant, GrayP^ = 55 ± 1 (Fig. 6), which is determined 
using the approach described in Section II C 2. 

It is noteworthy that the above calculations are based 
on the assumption that the dimer crystal of [32, 34] 
is the densest possible arrangement of hard tetrahedra. 
Although we cannot rule out the possibility that an 
even denser arrangement of tetrahedra that is different 
from the approximant and the dimer crystal might ex- 
ist, our observation that the dimer crystal is the dens- 
est structure that forms in simulations of 16 tetrahedra 
and fewer [32, 34] substantiates this assumption. The 
quasicrystal that we are using for comparison with the 
approximant has been assembled in simulations from the 
disordered fluid and therefore contains imperfections. We 
cannot rule out that a perfect quasicrystal might be ther- 
modynamically more stable than the approximant at all 



FIG. 6: Gibbs free energies of the approximant and the fluid 
close to the melting transition. The transition occurs at P^f — 
55. 

pressures. If this were the case, then the transition be- 
tween the approximant and the dimer crystal reported 
above would be substituted by a transition between the 
quasicrystal and the dimer crystal in the phase diagram. 
Therefore, while such a discovery could alter certain de- 
tails of the phase transition, it will not eliminate the ex- 
istence of a solid-solid phase transition reported in this 
work. 



D. Dimer-quasicrystal transformation 

To compare the relative thermodynamic stability of 
the dimer crystal and the quasicrystal in simulation, we 
set up a Monte Carlo simulation of a large dimer crys- 
tal with 2,916(= 4x9x9x9) tetrahedra in the iso- 
choric ensemble. To facilitate the transformation, the 
box dimensions are occasionally distorted in a random 
direction with the constraint that the total volume re- 
mains unchanged (variable-shape ensemble, [18]). This 
distortion allows the system to adjust to arbitrary lattice 
symmetries by relaxing shear stresses. 

We choose a constant packing density of <p = 50%, 
because at this density the quasicrystal is routinely ob- 
served to form spontaneously from the fluid. Structural 
changes are detected by counting the number of particles 
that form PDs and icosahedra using a shape-matching 
algorithm [47]; icosahdral motifs vanish when the qua- 
sicrystal forms [29]. Additionally, the pressure is deter- 
mined from the acceptance probability of trial volume 
changes as described in Section II B [17, 48]. 

The pressure shows a sharp spike after 4 million Monte 
Carlo cycles accompanying the melting of the dimer crys- 
tal (Fig. 7a). The spike quickly decreases to a plateau, 
which, after 15 — 20 million Monte Carlo cycles, relaxes 
to its equilibrium value. PDs and icosahedra form as the 
preferred local configurations in the melt (Fig. 7b). On 
the other hand, in the final solid structure, most particles 
are members of PDs and virtually no icosahedra remain. 
Diffraction images in Figs. 7c-f show that the final solid 
structure is the dodecagonal quasicrystal. The fact that 
the quasicrystal forms in the simulation with the melt as 
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FIG. 7: Transformation of the dimer crystal to the dodecago- 
nal quasicrystal in an isochoric simulation, (a) The pressure 
first spikes after 4 million Monte Carlo cycles and then relaxes 
during the melting of the dimer crystal. Between 15 and 20 
million Monte Carlo cycles, the quasicrystal forms from the 
melt, (b) The number of particles arranged in pentagonal 
dipyramids (PDs) or icosahedra (ico) increases rapidly dur- 
ing melting. In the quasicrystal essentially all particles form 
PDs while icosahedra disappear. Diffraction patterns confirm 
the transformation from the dimer crystal (c) to the melt (d,e) 
and then to the quasicrystal (f). 



an intermediate state confirms that both the quasicrys- 
tal and the melt are thermodynamically favored over the 
dimer crystal at the packing density <p = 50%. 
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FIG. 8: Relative stability of the dimer crystal, quasicrystal, 
and quasicrystal approximant. (a) Up to packing density 83% 
the dimer crystal has lower average free volume per particle. 
This helps to stabilize the approximant entropically. At high 
packing densities the dimer crystal should eventually have the 
highest average free volume, because its maximally achievable 
density is the highest of the three candidate structures, (b) 
Comparison of the Gibbs free energy differences between the 
dimer crystal and the approximant using the exact Frenkel- 
Ladd method, the mean-field approximation, and the cell- 
model approximations. The transition is predicted with all 
three methods even though the critical densities (j>A (approx- 
imant) and (j)D (dimer crystal) vary slightly. 



E. Origin of stability of the approximant 

To investigate the superior stability of the quasicrystal 
approximant compared to the dimer crystal over such a 
wide range of densities, we investigate the significance of 
collective particle motions by comparing the free energy 
estimates obtained from the mean-field approximation 
introduced in Section II D 3 with the exact free energy 
differences. We also analyze the dynamics in the ap- 
proximant by calculating the van Hove correlation func- 
tion [49] and visually inspecting the high-mobility parti- 
cles [50] in our simulations. 



1. Free volumes 

We calculate the mean-log average of the free volumes 
v f,ML of tetrahedra (Eq. 15) in the approximant, the 
dimer crystal, and the quasicrystal using the shooting 
method described in Section II D 1. The results are pre- 
sented in Fig. 8a. Whereas particles in the quasicrystal 
generally have a smaller v/,ml than in the approximant, 



the curves are shifted along the abscissa relative to one 
another by a fixed amount as indicated with arrows in 
Fig. 8a. This implies an identical thermodynamics for 
the quasicrystal and the approximant except for their 
different maximum packing densities. Indeed, tetrahedra 
experience similar local environments in the quasicrystal 
and its approximant. 

In contrast, the mean log free volume of the dimer 
crystal decays much more slowly with packing density 
and intersects the two other free volume curves. This 
finding suggests that the approximant relaxes more ef- 
ficiently during expansion, creating free volume for the 
particles more readily. Note that the packing density 
where the two curves cross is considerably below 84%, 
the density where the approximant becomes thermody- 
namically unstable, which underscores the significance of 
collective motions of particles in stabilizing the approxi- 
mant even at very high densities. 

The importance of collective motions may be further 
inferred by comparing the free energy difference esti- 
mated from a mean-field approximation with the exact 
value. As shown in Fig. 8b, the mean-field approximation 
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underestimates the stability of the approximant, which 
indicates that entropic contributions from collective mo- 
tions are significant. We suspect that slight rearrange- 
ments of particles in the approximant during expansion 
also increase its stability at lower packing densities. This 
is confirmed by estimating AG from a cell model ap- 
proximation. The cell model is similar to the mean-field 
approximation except that free volumes are calculated 
for a non-equilibrated structure obtained by isotropically 
expanding the densest packing to a given packing den- 
sity [45]. As shown in Fig. 8b, the cell-model approxima- 
tion underestimates the stability of the approximant even 
more than the mean-field approximation, which suggests 
the significance of small local rearrangements that occur 
while the structure is equilibrated after expansion. 



2. Dynamics in the approximant 

Correlated motions of tetrahedra are observed in long 
simulations of both the approximant and the quasicrys- 
tal at all densities. These motions are most apparent at 
packing densities below 65% where they give rise to lo- 
cal structural rearrangements, but they are still present 
at higher densities in the form of correlated vibrations 
of clusters of tetrahedra. The fundamental mechanism 
through which these rearrangements proceed is the rota- 
tion of single PDs around their principal axes by multi- 
ples of 72°. The rounded, disk-like shape of PDs, com- 
pared to tetrahedra with their sharp corners allows an 
easy rotation even in relatively dense configurations. 

The rotation of PDs is confirmed by observing several 
peaks in G s (r, t), the self-part of the van Hove correla- 
tion function [49], which implies that the tetrahedra in- 
deed move between discrete sites separated by geometric 
barriers (Fig. 9a, b). As reported in our earlier work [29], 
each tetrahedron in the quasicrystal and the approximant 
is part of a spanning network of interpenetrating PDs 
(that is, PDs that share a tetrahedron). The locations of 
the peaks in G s (r, t) correspond to the characteristic dis- 
tances of the nearest neighbor distances in the spanning 
network. 

We observe that not all PDs are equally likely to ro- 
tate. At high densities, the PDs capping the 12-fold rings 
in the center of logs (shown in green in Fig. lb) [29] rotate 
more frequently as they are spatially separated from the 
rest of the structure. This can be seen in the trajecto- 
ries of the high-mobility particles in the approximant at 
qb = 65% (Fig. 9c). Close to melting, however, rotations 
involve the full network of neighboring PDs, which allows 
the particles to diffuse over arbitrary distances (Fig. 9d). 
The underlying dynamics is identical in the quasicrystal. 
However the presence of defects leads to higher mobility 
in the quasicrystals that form in simulation as compared 
with "perfect" quasicrystals. Both the quasicrystal and 
the approximant exhibit some 'liquid-like' behavior since 
unlike simple crystals, diffusion can take place in these 
systems even in the absence of defects. 




r / a If 

FIG. 9: Particle dynamics in the quasicrystal approximant. 
(a,b) The self-parts of the van Hove correlation functions at 
4> = 65% (a) and = 50% (b) show various peaks, which indi- 
cates that the particles do not move continuously but have to 
overcome (geometric) barriers. The peak positions correspond 
to different levels of nearest neighbor distances in the under- 
lying PD network. (c,d) The trajectories of particles with the 
highest mobility are plotted. At high density, cj) — 65% (c), 
tetrahedra move along the edges of pentagons. This motion 
corresponds to rotations of the PDs in log centers. At in- 
termediate densities, = 50% (d), neighboring PDs start to 
rotate and the tetrahedra are more mobile. In the infinite 
time limit the tetrahedra can diffuse through successive PD 
rotations. 



At packing densities beyond 65%, PD rotations become 
extremely unlikely, but clusters of tetrahedra, including 
PDs, can still vibrate collectively. Figs. lOa-b show such 
correlated motions occurring in a time period of 50 mil- 
lion Monte Carlo cycles in a layer of the approximant 
at (j) = 75% and <p = 80% respectively. The vibrations 
are extremely slow, but their existence adds additional 
entropy to the system making the mean-field approxi- 
mation and the cell model inaccurate. No dynamics is 
observed in the dimer crystal. 

In general, thermodynamically equivalent local rear- 
rangements are a characteristic feature of quasicrystals 
and their approximant s. The transformation among 
these takes place via phason modes [51-53]. Elemen- 
tary excitations are phason flips, which previously have 
been observed with high-resolution transmission electron 
microscopy [54] and in simulations of two dimensional 
model systems [55]. 



IV. DISCUSSION AND CONCLUSION 

In general one might expect a 'simple' structure like the 
dimer crystal to form more easily than 'complex' struc- 
tures like the quasicrystal or its approximant. The ob- 
servation that tetrahedra defy this expectation suggests 
that structural complexity is not always a good indica- 
tor of thermodynamic stability. Indeed, although it has 
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FIG. 10: Correlated motion of clusters in a slab of the ap- 
proximant at (a) <j) — 75% and (b) <j) — 80%. Dark arrows 
correspond to the direction towards which each particle has 
moved after t = 5 X 10 7 Monte Carlo cycles; the length of 
each arrow is twice the distance the corresponding particle 
had travelled. There are several clusters of neighboring tetra- 
hedra moving collectively. A few of these clusters are high- 
lighted in blue. Not surprisingly, the mobility is higher at 
(j) — 75% as evidenced by longer arrows. 



been argued in the literature [35] that the dimer crys- 
tal first reported by [32, 34] and studied here might be 
the stable phase even at densities where the quasicrystal 
is reproducibly observed (down to densities of 50%), our 
free energy calculations demonstrate that the dimer crys- 
tal is in fact preferred thermodynamically only at very 
high densities (above 84%). On the other hand, inso- 
far as structural complexity increases a system's entropy, 
structurally complex arrangements of hard particles may 
be thermodynamically preferred over simpler ones. 

Indeed, we have shown that the structural features 
of the quasicrystal and the approximant allow for more 
complex dynamics than the dimer crystal at moderate 
and high densities as manifested in the behavior of the 
free volume as a function of packing density and the col- 
lective motions in the form of PD rotations. The exis- 
tence of the PD network facilitates collective particle mo- 
tions at low densities. Although rearrangements become 
vanishingly unlikely at higher densities, they appear to 
contribute additional entropy to the system and stabilize 
it over the dimer crystal, in which each particle can only 
'rattle' independently in its own cage. Rearrangements 
are impossible in the dimer crystal because no rearrange- 
able network exists there. 

The superior stability of the quasicrystal and its ap- 
proximant relative to the dimer crystal may also be at- 
tributed to the presence of almost-perfect face-to-face 
contacts between tetrahedra. There is a natural tendency 
for hard polyhedra to optimize face-to-face contacts at 
high densities in order to maximize configurational en- 
tropy. For instance, there are an infinite number of cubic 
arrangements of hard cubes with packing fraction one, 
but among them the simple cubic lattice, where all cubes 
are perfectly face-to-face, has the highest entropy and is 
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FIG. 11: Schematic phase diagram of hard tetrahedra sum- 
marizing our findings. In thermodynamic equilibrium the 
Dimer III crystal and the approximant are stable (Middle 
panel). In compression simulations the approximant is never 
observed, and only the quasicrystal forms. If crystallization 
is suppressed, then a jammed packing with local tetrahedral 
order forms [29, 36] (Lower panel). The transformation of the 
approximant or quasicrystal directly to and from the Dimer 
III crystal is not observed in simulation. Instead, during ex- 
pansion the Dimer III crystal transforms into the Dimer II 
crystal, and then the Dimer I phase prior to melting to the 
fluid (Upper Panel). 



thermodynamically stable [56]. 

Within the approximant, we observe that face-to- face 
contacts between neighboring tetrahedra are nearly per- 
fect in the sense that the touching faces are not signifi- 
cantly shifted with respect to one another. This is not 
true in the dimer crystal where inter-dimer face-to-face 
contacts are shifted and therefore not close to being per- 
fect. Abundance of strong face-to-face contacts makes 
the PD network more rearrangeable and collective mo- 
tions of particles more feasible, which in turn leads to a 
higher entropy and superior stability. 

We summarize our findings in a schematic phase dia- 
gram in Fig. 11. We note that hard tetrahedra are one of 
the few examples of hard particles with two distinct solid 
phases not mutually related by symmetry breaking. Our 
results show that entropic effects alone are sufficient for 
inducing highly nontrivial solid-solid phase transitions. 

Not all phase transformations are accessible in simu- 
lations on finite time scales. The observation that simu- 
lations only form the quasicrystal but never the approx- 
imant suggests that the quasicrystal is kinetically more 
easily accessible than the approximant - independent of 
whether it is thermodynamically preferred or not. This 
can be attributed to the fact that the transformation of 
a dodecagonal quasicrystal to one of its approximants 
proceeds through a process called zipper motion [57], 
which is extremely slow even in experiment [58]. Fur- 
thermore, transformation to the dimer crystal at packing 
densities greater than 84% is not observable in simula- 
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tions, and may be unobservable in experiments, due to 
the extremely slow kinetics at such high densities. 

In conclusion, we have shown that the quasicrystal 
and its approximant are thermodynamically favored over 
the dimer crystal at all experimentally realizable packing 
densities. We also observe a very rich dynamical behavior 
in the quasicrystal and its approximant induced by rota- 
tions of pentagonal dipyramids within an interconnected 
network. We have shown the significance of collective 
motions in stabilizing the approximant for a wide range 
of packing densities. 
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